#include "AnaDef.h"

TFile *ffile[]={0x0,0x0};
const Double_t EPSILON = 1e-12;

TString ffins[2];

//https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies#Tables_with_centrality_bins_for
const Double_t fncoll[]={1685, 1316,  921.2 ,  438.4 ,  127.7 ,  26.71};
//const Double_t fncsta[]={1.4e+02, 1.1e+02, 1.4e+02, 1.5e+02, 59, 18 , 4.4};
//const Double_t 

void getDoubleRatio(const TString tag)
{
  TList * ll = new TList; 

  const Int_t npar = 4;
  for(Int_t itype=0; itype<npar; itype++){
    TGraphAsymmErrors * gtmpStat[]={0x0,0x0};
    TGraphAsymmErrors * gtmpSys[]={0x0,0x0};
    for(Int_t ich = 0; ich<2; ich++){
      gtmpStat[ich] = (TGraphAsymmErrors*) ffile[ich]->Get(Form("grstat%d", itype));
      if(!gtmpStat[ich]){
        printf("no g ich %d itype %d\n", ich, itype); ffile[ich]->ls();exit(1);
      }
      gtmpSys[ich] = (TGraphAsymmErrors*) ffile[ich]->Get(Form("grsys%d", itype));
      if(!gtmpSys[ich]){
        printf("no g ich %d itype %d\n", ich, itype); ffile[ich]->ls();exit(1);
      }
    }

    TGraphAsymmErrors * gpnstat = new TGraphAsymmErrors;
    gpnstat->SetName(Form("grstat%d", itype));
    ll->Add(gpnstat);

    TGraphAsymmErrors * gpnsys = new TGraphAsymmErrors;
    gpnsys->SetName(Form("grsys%d", itype));
    ll->Add(gpnsys);

    /*
    Int_t np = gtmpStat[0]->GetN();
    if(np != gtmpStat[1]->GetN()){
      printf("check! not equal n! %d %d %s %s\n", np, gtmpStat[1]->GetN(), ffins[0].Data(), ffins[1].Data());

      const Double_t ch0x0 = gtmpStat[0]->GetX()[0];
      const Double_t ch1x0 = gtmpStat[1]->GetX()[0];
      if(fabs(ch0x0-ch1x0)>1e-4){
        for(Int_t itmpch=0; itmpch<2; itmpch++){
          for(Int_t itmpp=0; itmpp<gtmpStat[itmpch]->GetN(); itmpp++){
            const Double_t xx = gtmpStat[itmpch]->GetX()[itmpp];
            printf("please check ch %d ix %d xx %e\n", itmpch, itmpp, xx);
          }
        }
        exit(1);
      }
      else {
        if(gtmpStat[1]->GetN()<np){
          np = gtmpStat[1]->GetN();
          printf("check! correct np to %d\n", np);
        }
      }
    }
    */

    Int_t np = TMath::Min(gtmpStat[0]->GetN(), gtmpStat[1]->GetN());

    Double_t taai = 1;
    if(tag.Contains("RAA")){
      printf("setting taai!!");
      if(ffins[0].Contains("CEN0"))
        taai=fncoll[0];
      else if(ffins[0].Contains("CEN1"))
        taai=fncoll[1];
      else if(ffins[0].Contains("CEN2"))
        taai=fncoll[2];
      else if(ffins[0].Contains("CEN3"))
        taai=fncoll[3];
      else if(ffins[0].Contains("CEN4"))
        taai=fncoll[4];
      else if(ffins[0].Contains("CEN5"))
        taai=fncoll[5];
    }
    for(Int_t ip=0; ip<np; ip++){
      Double_t ix[]={-999,-999};
      Double_t iy[]={-999,-999};
      Double_t exl[]={-999,-999};
      Double_t exh[]={-999,-999};
      Double_t eyl[]={-999,-999};
      Double_t eyh[]={-999,-999};
      Double_t sysEyl[]={-999,-999};
      Double_t sysEyh[]={-999,-999};
      for(Int_t ich=0; ich<2; ich++){
        gtmpStat[ich]->GetPoint(ip, ix[ich], iy[ich]);

        exl[ich] = gtmpStat[ich]->GetErrorXlow(ip);
        exh[ich] = gtmpStat[ich]->GetErrorXhigh(ip);
        eyl[ich] = gtmpStat[ich]->GetErrorYlow(ip);
        eyh[ich] = gtmpStat[ich]->GetErrorYhigh(ip);
     
        sysEyl[ich] = gtmpSys[ich]->GetErrorYlow(ip);
        sysEyh[ich] = gtmpSys[ich]->GetErrorYhigh(ip);
      }

      if(fabs(ix[0]-ix[1])>2e-2){
        printf("not equal x! %d %e %e -- %e\n", ip, ix[0], ix[1], ix[0]-ix[1]); exit(1);
      }
      if(iy[0]<EPSILON || iy[1]<EPSILON){
        printf("negative y %e %e\n", iy[0], iy[1]); exit(1);
      }
      const Double_t dr = iy[0]/iy[1]/taai;
      const Double_t erryl = dr * TMath::Sqrt(TMath::Power(eyl[0]/iy[0],2) + TMath::Power(eyl[1]/iy[1],2) );
      const Double_t erryh = dr * TMath::Sqrt(TMath::Power(eyh[0]/iy[0],2) + TMath::Power(eyh[1]/iy[1],2) );
      
      const Double_t sysErryl = dr * TMath::Sqrt(TMath::Power(sysEyl[0]/iy[0],2) + TMath::Power(sysEyl[1]/iy[1],2) );
      const Double_t sysErryh = dr * TMath::Sqrt(TMath::Power(sysEyh[0]/iy[0],2) + TMath::Power(sysEyh[1]/iy[1],2) );
      
      printf("ip %d %e %e %e %e -- %e, %e %e %e %e -- %e %e %e %e -- taai %e\n", ip, ix[0], ix[0]-ix[1], iy[0], iy[1], dr, exl[0], eyl[0], eyl[1], erryl, sysEyl[0], sysEyl[1], sysErryl, sysErryh, taai);

      gpnstat->SetPoint(ip, ix[0], dr);
      gpnstat->SetPointError(ip, exl[0], exh[0], erryl, erryh);

      gpnsys->SetPoint(ip, ix[0], dr);
      //too large to see the points
      gpnsys->SetPointError(ip, exl[0], exh[0], sysErryl, sysErryh);
    }

  }

  TFile * fout=new TFile(tag+".root","recreate");
  ll->Write();
  fout->Save();
  fout->Close();
  delete fout;

}

int main(int argc, char * argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d %s\n", ii, argv[ii]);
  }

  if(argc!=3){
    printf("wrong argc!\n"); return 1;
  }

  TString pwd=gSystem->pwd();
  pwd.ToUpper();

  TString tag(pwd.Contains("RAA")?"RAA":"DR"); 

  for(Int_t ich=0; ich<2; ich++){
    TString fin(argv[ich+1]);
    ffins[ich]=argv[ich+1];

    ffile[ich] = new TFile(fin);
    if(!ffile[ich]){
      printf("no file %d %s\n", ich, fin.Data());exit(1);
    }

    tag+="_";
    fin=fin(fin.First('/')+1, fin.Length());
    fin=fin(0,fin.First('.'));
    tag += fin;
  }

  getDoubleRatio(tag);

  return 0;
}
